#include <stdio.h>
#include <stdlib.h>
#include "gmp.h"

//Parameters for GMP lib
#define ITER            25
#define PRECISION       100000000
#define BASE            10
#define DIGITS          30000000

//Calculates an1 =  ((a0 + b0) / 2.0);
void an1(mpf_t, mpf_t, mpf_t);

//Calculates bn1 = sqrt(a0 * b0);
void bn1(mpf_t, mpf_t, mpf_t);

//Calculates tn1 = tn - pn(pow(an-an1,2))
void tn1(mpf_t, mpf_t, mpf_t, mpf_t, mpf_t, mpf_t);

//Calculates  pn1 =  2.0 * p0
void pn1(mpf_t, mpf_t);

//Main function
int main(int argc, char *argv[]) {
//Set new iteration number
  int iter;
    if(argc == 2)
        iter = atoi(argv[1]);
    else iter = ITER;

    fprintf(stderr,"Calculando Gauss-Legende Sequencial com %d iteracoes...\n",iter);

    //Declaration of main variables
    mpf_t a0, b0, t0, p0, a1, b1, t1, p1;

    //Declaration of aux variables
    mpf_t aux, result, tmp;

    //Define precision
    mpf_set_default_prec(PRECISION);

    //initialization of variables
    mpf_inits(a0, b0, t0, p0, a1, b1, t1, p1, result, tmp, aux, NULL);

    //Set a0 initial valor
    //    a0 = 1.0;
    mpf_init_set_ui(a0, 1);

    //Set b0 initial valor
    //    b0 = 1.0 / M_SQRT2;
    mpf_init_set_ui(b0, 1);
    mpf_init_set_ui(tmp, 2);
    mpf_sqrt(tmp, tmp);
    mpf_div(b0, b0, tmp);

    //Set t0 initial valor
    //    t0 = 1.0 / 4.0;
    mpf_init_set_ui(t0, 1);
    mpf_init_set_ui(tmp, 4);
    mpf_div(t0, t0, tmp);


    //Set p0 initial valor
    //    p0 = 1.0;
    mpf_init_set_ui(p0, 1);


    //Iterations of a,b,t,p
    int i = 0;
    for (i; i < iter; i++) {

        //Calculus of each iteration of a,b,t,p
        an1(a0, b0, a1);
        bn1(a0, b0, b1);
        tn1(a0, a1, b0, p0, t0, t1);
        pn1(p0, p1);


        //Set new values to old variables to continue the iterations
        //        a0 = a1;
        mpf_set(a0, a1);
        //        b0 = b1;
        mpf_set(b0, b1);
        //        p0 = p1;
        mpf_set(p0, p1);
        //          t0 = t1;
        mpf_set(t0, t1);



    }

    //Calculus of pi using previous calculated variables 
    // pi = pow(an+bn,2)/(4*tn)
    mpf_add(tmp, a1, b1);
    mpf_pow_ui(tmp, tmp, 2);
    mpf_mul_ui(aux, t1, 4);
    mpf_div(result, tmp, aux);


    //Print the result
    mpf_out_str(stdout, BASE, DIGITS, result);
    //Clear variables
    mpf_clears(a0, b0, t0, p0, a1, b1, t1, p1, aux, result, tmp, NULL);


    return 0;
}


// an1 =  ((a0 + b0) / 2.0);
void an1(mpf_t a0, mpf_t b0, mpf_t an1) {

    mpf_add(an1, a0, b0);
    mpf_div_ui(an1, an1, 2);


    return;
}


//bn1 = sqrt(a0 * b0);
void bn1(mpf_t a0, mpf_t b0, mpf_t bn1) {

    mpf_mul(bn1, a0, b0);
    mpf_sqrt(bn1, bn1);


    return;


}


//tn1 = tn - pn(pow(an-an1,2))
void tn1(mpf_t a0, mpf_t a1, mpf_t b0, mpf_t p0, mpf_t t0, mpf_t tn1) {



    //   tmp = pow(a0 - a1, 2);
    mpf_sub(tn1, a0, a1);
    mpf_pow_ui(tn1, tn1, 2);

    //    (t0 - (p0 * (tmp)))
    mpf_mul(tn1, p0, tn1);
    mpf_sub(tn1, t0, tn1);


    return;
}


//  pn1 =  2.0 * p0
void pn1(mpf_t p0, mpf_t pn1) {

    //    2.0 * p0
    mpf_mul_ui(pn1, p0, 2);

    return;
}
